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Abstract 

Background: Analyzing and monitoring uterine contractions during pregnancy is relevant to the field of 
reproductive health assessment. Its clinical importance is grounded in the need to reliably predict the onset of labor 
at term and pre-term. Preterm births can cause health problems or even be fatal for the fetus. Currently, there are no 
objective methods for consistently predicting the onset of labor based on sensing of the mechanical or 
electrophysiological aspects of uterine contractions. Therefore, modeling uterine contractions could help to better 
interpret such measurements and to develop more accurate methods for predicting labor. In this work, we develop a 
multiscale forward electromagnetic model of myometrial contractions during pregnancy. In particular, we introduce a 
model of myometrial current source densities and compute its magnetic field and action potential at the abdominal 
surface, using Maxwell's equations and a four-compartment volume conductor geometry. To model the current 
source density at the myometrium we use a bidomain approach. We consider a modified version of the Fitzhugh- 
Nagumo (FHN) equation for modeling ionic currents in each myocyte, assuming a plateau-type transmembrane 
potential, and we incorporate the anisotropic nature of the uterus by designing conductivity-tensor fields. 

Results: We illustrate our modeling approach considering a spherical uterus and one pacemaker located in the 
fundus. We obtained a travelling transmembrane potential depolarizing from —56 mVto — 16 mVand an average 
potential in the plateau area of —25 mV with a duration, before hyperpolarization, of 35 s, which is a good 
approximation with respect to the average recorded transmembrane potentials at term reported in the technical 
literature. Similarly, the percentage of myometrial cells contracting as a function of time had the same symmetric 
properties and duration as the intrauterine pressure waveforms of a pregnant human myometrium at term. 

Conclusions: We introduced a multiscale modeling approach of uterine contractions which allows for incorporating 
electrophysiological and anatomical knowledge of the myometrium jointly. Our results are in good agreement with 
the values reported in the experimental technical literature, and these are potentially important as a tool for helping in 
the characterization of contractions and for predicting labor using magnetomyography (MMG) and 
electromyography (EMG). 



Background 

Modeling the myometrium contractility during pregnancy 
is of clinical importance, since it can aid in understand- 
ing the mechanism of labor, and, thus, it can help in 
monitoring the health of both the fetus and the mother. 
The occurrence of labor is in general accompanied by 
the appearance of periodic contractions which increase 
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the intrauterine pressure to the point that cervix dilata- 
tion is manifested [1]. However, from clinical experi- 
ences, not all uterine contractions are efficient, i.e., lead 
to labor. Term labor is expected to occur after the 37th 
week of pregnancy, but in the last decade the incidence 
of preterm labor has increased significantly (12% of all 
births). Preterm birth can cause health problems or even 
be fatal for the fetus if it happens too early, and, at the same 
time, it imposes significant financial burdens on health 
care systems [2]. Therefore, it becomes critical to better 
understand the mechanism of bioreproduction which, as 
a consequence, would allow for the development of more 
effective forms of therapy that might help to predict labor 
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and control the occurrence of labor. In the following we 
will discuss briefly the uterine microanatomy, previous 
contractions models, and our modeling approach. 

Uterine microanatomy 

The adult uterus is a thick walled, hollow, muscular organ 
formed by three layers: the external serous perimetrium, 
the myometrium, and the inner mucous endometrium [1]. 
The non-pregnant uterus wall thickness is approximately 
15 mm to 20 mm, and during pregnancy the uterine 
wall becomes thin with values at the 39th week of preg- 
nancy ranging from 7.4 ±1.8 mm at the low anterior 
wall (lower segment) at the bladder interface, to 10.06 
zb 1.9 mm at the posterior wall [3]. The myometrium 
is responsible for contractions, and it is formed by fas- 
ciculi which are comprised of sheet-like and cylindrical 
bundles of myocytes embedded in a connective tissue 
matrix [4] . The myocytes in a cylindrical bundle contract, 
thus shortening the smooth tissue and increasing uterus 
wall tension, hence increasing the intrauterine pressure. 
Figure 1 illustrates the microanatomy of the pregnant 
human myometrium. 

The uterine microanatomy is consistent with action 
potential propagation [4]: (i) myocytes are densely packed 
within a bundle, (ii) bundles are contiguous within a fasci- 
culus, and (iii) fasciculi are contiguous via communicating 
bridges formed with myocytes. In addition, the uterine 
changes during gestation are accompanied by the forma- 
tion of gap junctions, which are one of the mechanisms 
for transmitting of contractile activity from cell to cell in 
a coordinated manner [1,4]. The structure of the fasicu- 
lata within the uterus has not yet been well defined, but 
generally it makes the propagation of the action potential 
anisotropic [5,6]. 

Uterine contraction models 

Uterine contractions can be described by their mechan- 
ical and electrophysiological aspects. A mechanical 



contraction is manifested as a result of the excitation as 
well as the propagation of electrical activities in the uter- 
ine muscle, and appears in the form of an intrauterine 
pressure increase. Existing models approach the problem 
separately at the organ level [7-12] or at the cellular level 
[13-17]. 

At the organ level, the models presented in [7-11] focus 
on predicting the contractile forces that closely resem- 
ble clinical measurements of normal intrauterine pres- 
sure during contractions in labor, and, more recently, the 
action potential propagation. In [7], the authors assume 
that the uterus is a hollow ovoid formed by discrete 
contractile elements that propagate electrical impulses, 
generate tension, and have defined contracting and refrac- 
tory periods. The envisioned mechanism for intercellular 
communication is based on action potential propagation, 
which is simulated by using a discrete state model for each 
cell. In [8], the authors revisit the model developed in [7] 
and perform multiple experiments with different uterine 
shapes, cell numbers, and initial distributions of active 
and resting cells. In [9], the author uses a discrete state 
model for combining action potential propagation and 
intercellular calcium wave propagation, two mechanisms 
of intercellular communication. However, in [7-9], math- 
ematical and physical descriptions of the models are not 
provided. In [11], the trajectories of growth and myome- 
trial tension at the onset of labor are modelled using a 
statistical modeling approach. The authors model the tra- 
jectories of uterine wall thickness, volume, and tension 
as functions of gestational age using a prolate ellipsoid 
method, intrauterine pressure results from the literature, 
and ultrasound measurements of the shape of the uterus 
collected on 320 subjects. Regarding the electrical activ- 
ity of uterine contractions, a model of myometrial action 
potential propagation as measured by surface electrohys- 
terography is proposed in [10]. The authors develop a 
myometrium skin conduction model consisting of four 
parrallel layers (myometrium,abdomen,fat, and skin), and 
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model the action potential using a Gamma probability 
distribution. This model assumes that the electrical con- 
ductivity of the myometrium is isotropic and that the 
abdominal curvature is negligible, thus making it suitable 
to model the electrical potential in a limited area of the 
abdominal surface. Recently, in [12], the authors model 
the electrical propagation of action potential using a 3- 
dimensional myometrium for different initial conditions 
and stimuli. The authors use a reaction-difussion equation 
coupled with a Fizhugh-Nagumo model of ionic currents 
and consider a homogeneous isotropic myometrium. 

At the cellular level, the models focus on predicting 
the changes of ionic concentrations in the intracellular 
and extracellular media during a contraction, and, as a 
consequence, on modeling the transmembrane potential 
evolution of a myocyte as a function of time. In [13,14] 
a model is developed to simulate the complete process 
of a single myometrial smooth muscle contraction, which 
is initiated by depolarization. The model is based on the 
electrophysiological properties of a myocyte and on the 
cellular mechanisms that relate the rise in concentration 
levels of intracellular ion calcium Ca 2+ to stress produc- 
tion. In [17], the authors model all known individual ionic 
currents of myometrial smooth muscle close to labor and 
combined them into a mathematical model of myometrial 
action potential generation. The model is shown to suc- 
cessfully mimic several recordings of spontaneous AP and 
force in uterine smooth muscle. 

In this work, we propose a multiscale forward elec- 
tromagnetic model of human myometrial contractions 
during pregnancy that jointly takes into account electro- 
physiological and anatomical knowledge at the cellular, 
tissue, and organ levels. Our model aims to help in the 
characterization of contractions and the prediction of 
labor using MMG [18] and EMG [19]. Here, we extend 
our partial results presented in [20]. Figure 2 illustrates 
the different levels considered in our modeling approach. 
In particular, our approach is twofold: first, we model the 



current source density at the myometrium, using models 
of myocyte electrophysiological activity and anisotropic 
conductivity, and second, we solve the forward electro- 
magnetic problem, specifically, we compute the magnetic 
field and the action potential at the abdominal sur- 
face generated by the myometrial current-source density 
using Maxwells equations subject to a volume conduc- 
tor geometry. To model the current source density at the 
myometrium we propose to apply a bidomain approach. 
The bidomain equations are a set of reaction-diffusion 
equations derived first for modeling the current sources 
of the myocardium as a function of the cardiac-myocyte 
transmembrane potential, and these equations proved to 
be a successful approach to study the functioning of the 
heart [21,22]. The diffusion part of the equations gov- 
erns the spatial evolution of the transmembrane potential, 
and the reaction part is given by the local ionic current 
cell dynamics. Here we introduce a modified version of 
the FitzHugh-Nagumo (FHN) equation for modeling ionic 
currents in each myocyte. Though FHN does not con- 
sider explicitly the Ca 2+ dynamics, the simplicity of the 
FHN model makes it an attractive candidate for model- 
ing the propagation of depolarization waves in such large 
2D and 3D simulations as the numerical examples pre- 
sented in this work. We propose a general approach to 
design the conductivity tensor orientation for any uterine 
shape and estimate the conductivity tensor values in the 
extracellular and intracellular domains, using Archies law 
[23]. We illustrate our modeling approach using a spher- 
ical, four-compartment volume conductor geometry. As 
we elaborate later on, certain aspects of it stand as a fair 
approximation of the volume conductor geometry when 
performing MMG measurements. 

The notational convention adopted in this paper is as 
follows: italic font indicates a scalar quantity, as in a; low- 
ercase boldface indicates a vector quantity, as in a, except 
for vector fields used in Maxwells equations such as elec- 
tric field E, magnetic field B, and current density /; upper 
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case italic indicates a matrix quantity, as in A. The matrix 
transpose is indicated by a superscript " r " as in A 1 ', and 
the identity matrix of size n x n is denoted I n . The set S> n 
denotes the vector space of symmetric nxn matrices, and 
the spaces of nonnegative definite matrices and positive 
definite matrices are denoted by and respectively. 
The inner product and norm denned in the Euclidean 
space are denoted by (•, •) and ||-||, respectively. 

Methods 

In this section we discuss the electromagnetic source rep- 
resentation of uterine contractions. We will introduce a 
four-compartment volume conductor model formed by 
an anisotropic bidomain myometrium and will present 
the models for the extrauterine electrical potential, mag- 
netic field, and myometrial current source density. Finally, 
we will present a numerical example to illustrate our 
modeling approach. 

Volume conductor model 

Figure 3 illustrates the four-compartment volume conduc- 
tor geometry for our problem, where A represents the 
abdominal cavity and dA the boundary surface denned by 
the abdomen, M represents the myometrium, and dM 
and dU are its external and internal boundary surfaces, 
respectively. The volume denoted by U represents the 
space filled with the amniotic fluid that exists between the 
internal uterine wall dU and the boundary dT denned by 




Figure 3 Representation of the four-compartment volume 
conductor geometry and the forward electromagnetic problem 
of uterine contractions. 



the fetus volume T . The vectors r and r indicate the posi- 
tions of the observation point and source, respectively, 
with respect to the main axis of reference. 

Extrauterine magnetic field and electrical potential 

The electromagnetic properties of uterine contractions 
can be analized by solving a set of Maxwells equations 
[24,25] subject to boundary conditions given by the vol- 
ume conductor geometry. In general, in a passive non- 
magnetic medium the total current density is the sum 
of the ohmic volume current and the polarization cur- 
rent, also known as the displacement current [24,25]. The 
ohmic current is the result of ions flowing in the medium 
and the displacement current is the result of a time vary- 
ing electric field [24,25]. If the temporal variations of the 
electric field and the magnetic field, i.e. low frequencies, 
are small enough that the displacement current is very 
small with respect to the ohmic current, then it is valid to 
model the electromagnetic phenomena using the quasi- 
static approximation of the Maxwell Equations. Under the 
quasi-static approximation, changes in the sources that 
generate the electromagnetic field affect all field quanti- 
ties instantaneously in the whole domain, and the total 
current density in the volume conductor geometry will be 
the sum of ohmic currents only. Biolectromagnetic fields 
vary slowly in time, with frequency components below 
lKHz and with a spatial characteristic length scale several 
times much larger than the volume conductor of inter- 
est [26], i.e, much larger than the diameter of the uterus. 
Hence, changes in the bioelectric sources affect the bio- 
electromagnetic field in the volume conductor geometry 
instantaneously, which justifies the use of the quasi-static 
approximation of Maxwells equations [25,26]. Therefore, 
the extrauterine magnetic field B(r, t) at a position r and 
instant t is given as follows: 

Vx —B(r,t)=J(r 1 t), (1) 

Mo 

where fi Q is the permeability of free space and / (r, t) is the 
total current density (in A/m 2 ). / (r, t) is given by 

J(r,t)=J s (r,t) + G(r)E(r,t), (2) 

where / s (r, t) is the uterine current density source and 
G (r) E (r, t) is the conduction current density (or return 
currents), as described by Ohms law, E(r>t) denoting 
the electric field established by / s (r, t) and G (r) e S^_ + 
denoting the conductivity tensor denned by each com- 
partment. Then, from the quasi-static conditions, V • 
/ ( r , t) = 0, so V • G (r) E(r,t) = - V / s (r, t). Moreover, 
since V x E (r, t) = 0, it follows that E (r, t) = - V0 (r, t), 
where (p (r, t) denotes the potential. Thus, the equation 
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that governs the relationship between electromyogram 
potentials and uterine current sources is 

V-G(r)V0(r,O = V-/ 8 (r,0. (3) 

Solving the forward electromagnetic problem of uter- 
ine contractions requires the computation of B(r, t) and 
(j) (r, t) at dA using Eqs.(l) and (3), assuming that / s (r, t) 
is known in M and G (r) is known in all domains denned 
by the volume conductor geometry (see Figure 3). 

The biological current sources / s (r, t) in the myo- 
metrium are the transmembrane ionic fluxes due to 
concentration gradients, which flow across the surface 
membrane of the myocyte (smooth cells) from the extra- 
cellular medium into the intracellular medium and vice 
versa. The density of these ionic currents is also referred 
to as the impressed current density, since its origin is non 
electrical in nature, and it is the primary cause for the 
establishment of an electric field that induces secondary 
density currents in a conductive domain. We will model 
/ s (r, t) using a bidomain approach, which has proved to 
be a successful method to study electrophysiological activ- 
ity in the myocardium [21,22] and, more recently, in the 
uterus [27]. 

Myometrial current-source density model 

In the myometrium, the intracellular and extracellular 
domains are both physically connected through mem- 
brane gates, and the intracellular domain is connected 
though gap junctions [1,19]. Therefore, we model the 
myometrium using the bidomain modeling approach. 
This approach represents the tissue (myometrium) as two 
interpenetrating extra-intracellular continuous domains 
with different conductivity values along and across the 
direction of the fiber [21,22], and it models the tissue 



using the generalized-passive cable equation. The bido- 
main modeling approach was originally derived for mod- 
eling the propagation of the transmembrane potential of 
the myocardium and proved to be a successful approach 
to study the functions of the heart [21,22]. Figure 4 shows 
a simplified illustration of the tissue and the bidomain 
approach, where <p[ (r, t) and 0 e (r> t) are the intracellu- 
lar and interstitial potentials, respectively, and v m (r, t) = 
0i (r, t) — 0 e (r, t) is transmembrane potential. The con- 
ductivity tensors in the intracellular and extracellular 
domains are denoted by G { and G e (in S/m), and, using 
Ohms law, the current densities in each domain are given 
by /i >e (r, t) = —G [e V0i >e (r, t) . The transmembrane vol- 
ume current density in (A/m 3 ) is denoted by j m (r, t) and 
is given by 



jm {f>t) — Cl m C m + /i on /stim> (4) 

at 

( 9Vm \ , . 

— a m I c m r. ~r/ion /stim I > WJ 



where j[ on (r, t) is the ionic volume current density (in 
A/m 3 ) of a myocyte, y'stim (r, t) is the stimulus volume cur- 
rent density (in A/m 3 ), c m is the membrane capacitance 
per unit area (in F/m 2 ), and a m is the surface-to-volume 
ratio of the membrane (in 1/m). To summarize the above 
description, the bidomain approach models the signal 
propagation in the tissue using the generalized passive 
cable equation model (Figure 4) which considers the cur- 
rent densities flowing through the extracellular domain 
given by / e (r, t), flowing through the intracellular domain 
given by J x (r, t), and connecting both media given by 
jm (r, t). 




Figure 4 Illustration of the bidomain modeling approach. 
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Applying the conservation of charge to both domains, 
we obtain the following relationships: 

V'h(r,t)=j m (r,t), and (6) 
V-Ji(r,t) = -j m (r,f). (7) 

Adding (6) and (7), we have V • (J { (r, t) + / e (r, t)) = 

0. Hence, the total current density in the myometrium is 
given by 

/ (r, t) = -G[ V0i (r, t) - G e V(j) e (r, t) , r e M, (8) 

which can be expressed in terms of v m (r, t) and (j) e (r, t) as 
follows: 

/ (r, t) = -G[ Vv m (r, 0 - G M V0 e (r, t) , r e M, (9) 

where G M = (G e + Gj^ e §+ + is the bulk myometrium 
conductivity tensor. Since spatial variations of v m (r, t) 
depend on the local establishment of a transmembrane 
current density, y m (r, t) ^ 0, we define the impressed 
current-density source as / s (r, t) = — G|Vv m (r, t). Note 
that / s (r, £) exists only when the spatial gradient exists, 

1. e., only in a region where the myometrium is undergoing 
depolarization (excitation) or repolarization. 

The total current at the myometrium / (r, t) depends 
on the spatio-temporal variations of v m (r, t) and <p e (r, t), 
which are governed by the system of equations formed by 
Eqs. (5), (6) and (7). Using simple algebraic manipulations, 
the aforementioned system of equations can be written in 
terms of v m (r, t) and 0 e ( r > t) only, obtaining the following 
equivalent expressions: 

/ / 9v m (r, t) 

V • GjV(v m (r, t) + 0 e (r, t)) = a m lc m m ^ 

+/ion(^ 0 -/stimO, ^ 

(10) 

V • (G[ + g;> V0 e (r, t) = - V • G(Vv m (r, t) . (11) 

This set of reaction-diffusion equations is also known 
as the bidomain equations [21,22]. The diffusion part of 
the equations governs the spatial evolution of both the 
transmembrane and extracellular potentials, and the reac- 
tion part is given by the local ionic current cell dynam- 
ics. The solutions for v m (r, t) and 0 e (r> t) depend on 
/ion ( r > 0>/stim ( r > t)> an d the conductivity tensors, in addi- 
tion to the boundary and the initial conditions. Since 
our goal is to model the propagation of electrical activ- 
ity in the myometrium, we are interested in the class of 
traveling-wave solutions of these equations whose wave- 
form depends on /i on (r, t) and whose initiation depends 
on /stim (r, t). In what follows, we describe the models for 
both current densities, J lon (r, t) and / s ti m ( r > t)> an d for the 
conductivity tensors, G[ and G^. 



Ionic current model 

The predominant type of transmembrane-potential wave- 
forms measured in the pregnant human myometrium 
are spikes and plateau [1,28-31]. In this work, we focus 
on modeling the plateau-type transmembrane potential, 
as it has been more frequently observed [1,28-30]. In 
particular, we model J lon (r, t), using a variation of the 
FitzHugh-Nagumo (FHN) equations [32-34], as follows: 

/ionO> t) = ~— (I< (V m - Vi) (V2 - V m ) 

x (v m - v 3 ) - w) , and (12) 

dw 

— =€ 2 Wv m -yw + 8), (13) 
dt 

where ei, 62, k, v\, V2, V3, 5, y, and ft are model constants, 
and w (in V) is a state variable of the model. The parameter 
61 (in £2m 2 ) controls the sharpness of the leading and trail- 
ing edges of the action potential waveform; the smaller e\ 
is, the more vertical the edge is. Note that e\ is a quantity 
of resistivity, therefore the smaller its value, the greater the 
permeability of the membrane to ionic flux. The parame- 
ter 62 (in s _1 ) controls the action potential duration; the 
smaller 62 is, the longer it takes a cell to recover. The 
parameters vi, V2, V3 (in V), and k (in 1/V 2 ) control the 
range of v m (r, t). Note that for a given set of /<, vi, V2, and 
V3, the ratios ft/y and 8/y control the excitability thresh- 
old of the cell. The larger /3/y is, the lower the excitability 
threshold that sets the cell dynamic to an oscillatory sta- 
ble behavior between resting and exciting states is. Over 
a certain value, the cell dynamic becomes bistable, that 
is, if the cell starts from a resting potential, it changes to 
an excited state and remains there. On the other hand, 
a very negative /3/y value results in a permanent rest- 
ing state. In the Results Section, presented further below, 
we select the model parameters using phase-space anal- 
ysis and using the transmembrane potentials recorded 
from isolated human myometrial strips at term as a refer- 
ence [13,30]. This model does not consider explicitly the 
Ca 2+ dynamics, and, moreover, it assumes that changes 
in the intra- and extra-cellular ion concentrations are 
insignificant even after several depolarizations. However, 
its simplicity facilitates the modeling of the propagation of 
depolarization waves in large 2D and 3D domains. 

Stimulus current model 

We also introduce a temporal-spatial model for J st [ m > rep- 
resenting the stimulus due to pacemaker areas [1,19], as 
follows: 




La Rosa etal. BMC Medical Physics 2012, 12:4 
http://www.biomedcentral.eom/1 756-6649/1 2/4 



Page 7 of 16 



where /z;(r, t) is a spatio-temporal function with range in 
[0, 1], Vi is the amplitude (in V), and N p is the number of 
pacemaker areas. Intuitively, the former should modify the 
excitability of the cell at a certain instant of time based on 
the threshold value. In particular, our model assumes that 
the uterine myocyte can act as either a pacemaker or pace- 
follower, specifically, the spontaneous electrical behavior 
exhibited by the myometrium is an inherent property of 
the uterine myocyte (see [19] for more details.) Note that 
the size, duration, and intensity of the pacemaker area 
need to be chosen such that a stable traveling waveform 
solution to the bidomain equations on the myometrium is 
possible. 

Conductivity tensor model 

The structure of the fasiculata within the uterus has not 
been completely characterized as a whole. Despite of not 
having a sense of a global uterine-fiber arcuitecture, it has 
been possible to characterize local structures in each of 
the three layers of the uterine wall [5,6,12]. In [5], the 
authors investigated the global fiber architecture of the 
non-pregnant uterus by diffusion tensor magnetic reso- 
nance imaging (DTMRI). From the ex- vivo analysis of five 
non-pregnant uteri, the authors identified an inner circu- 
lar layer around the uterine cavity on slices orthogonal 
to the long axis of the organ. In the regions outside the 
inner circular layer, they could not identify a global struc- 
ture but did find several locally aligned groups of fibers. At 
the level of the cervix, they found an outer circular layer 
and an inner region with mostly longitudinal components. 
In [12], from the analysis of the DTMRI of one postpar- 
tum uterus, the authors concluded that the myometrium 
is organized into bundles, with the fibre bundles forming 
interweaving sheets. In the following, we will introduce 
an approach for designing the conductivity tensors in the 
myometrium. 

Assume that the conductivity tensors are diagonal in 
a local coordinate system that is defined with respect 
to each myocyte and characterized by the unit vectors 
{ei,C2> e 3}- in particular, G[ and G e are diagonal matrices 
e §+ + given by 



/cr ix 0 0 \ 

0 G [y 0 

V 0 0 cr iz J 



/a ex 0 0 \ 
0 0r ey 0 

V 0 0 o e J 



In order to take into account variable fiber orientation in 
the myometrium, we need to describe it in a global Carte- 
sian coordinate system in which the local basis is defined 
at any point r as A= [a\(r) a<i(r) 03(1*)], where a<$(r) is 
parallel to the main fiber axis. The representation of the 



tensors G\ and G e in terms of a global coordinate system 
is given by 



G[ = AG[A T , G ! e =AG e A T . 



(16) 



Assuming that the myocyte fiber conductivities in both 
domains are cylindrically symmetric, then cr ex = o ey = 
a et , <7i x = <7i y = ait, <7[ z = an, and <r ez = <r e \. Therefore, the 
conductivity tensors can be expressed as follows [35] : 

G[ = fel - tfit) a 3 (r)a^(r) + cr it I 3 , and (17) 
G f e = (orei - cr et ) a 3 (r)a%(r) + cr et / 3 . (18) 

Hence, to construct the conductivity tensors as a func- 
tion of r, it is enough to define the vector field a^r), 
the myometrial fiber orientation, in each location of the 
anisotropic domain, as well as the conductivity values 
(Jii, cr e i, errand a et , because of the assumption of cylindric 
symmetry. 



Designing myometrial fiber orientations To design 
a^if) at each point r, we represent the uterus as a hol- 
low volume with uniform thickness, and we describe it as 
a union of mutually disjoint closed surfaces or layers. We 
use the implicit definition of a surface, namely, the set of 
points r satisfying /(r) = 0. Then, at each point r, we 
define a set of local orthonormal coordinate axes given 
by \n(r),t\(r), t2(r)}> where n(r) = ||v/(r)|| * s tne nor " 
mal vector to the layer containing r, and t\(r) and t^if) 
are mutually orthogonal vectors that belong to the tangent 
plane of the respective layer at point r. We define t\ (r) and 
t2(r), using the curve of symmetry of the uterine inner- 
circular layer as a reference [5] . This curve goes from the 
fundus to the cervix, and it coincides with the long axis 
of the non-pregnant uterus. We denote C as the curve of 
symmetry using the following parametric representation 
as a function of a single parameter t: 



C:t\-^ r c (t), t 1 <t<t 2 , 



(19) 



where yq it) is a point defined with respect to the 
global coordinate system, and yq (t\) and yq (£2) are the 
extreme points of the curve. For example, let yq it) = 
(x it) , y it) , z it)) with respect to the Cartesian system. 



(15) Define k(r) = 



dr c {t) 
dt 



dr c {t) 
dt 



the unitary vector 



field with direction given by the tangent vector of yq it) 
at to, where yq (£0) is the closest point to r such that 

/ dr c {t) 
\ dt 



> ?c ito) r) = 0. Then, we define t\ (r) to be con- 



tained in the plane formed by /<(r) and n(r) as follows: 



fi(r) = pk(r) + yn(r), 



(20) 
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subject to the following conditions: 
(?i(r),«(r)) = 0, and 



(?i(r),?i(r)) 2 



1. 



(21) 
(22) 

Therefore, replacing (20) in (21) and (22) we obtain the 
following system of equations: 

P |£(r),/z(r)j + y = 0, and (23) 

P 2 + y 2 + 2 0 y (*(r), *(r)) = 1. (24) 

Solving for all r, such that ^(r),w(r)j 7^ 1, we obtain 
P = ±- 



(£(r),»(r)) 

and y = =F — — 



= . Then £2 W = 

(?(r),*(r)) " ' /l-(?(r),«(r)) 

£i (r) x n(f) = /3k(r) x n(r), since by definition it is mutu- 
ally orthogonal to t\(r). Hence, given t\(r) and ti(y)> we 
define a<$(r) as follows: 



a^ir) = t\(r) cos (a) + £2(7*) sm (°0 » 



(25) 



where a is the fiber orientation angle with respect to £1 (r). 
In order to take into account complex fiber orientations 
a can be modeled as a spatial function denned over the 
domain of interest. Given our uterine volume assump- 
tions, the points tq (£1) and tq (£2) are the only points that 
satisfy the condition |/c(r), #f(r)j = 1. Since at rc (h) and 

*"c (£2) we cannot define t\(r) and £2(7*) using the curve C 
as a global reference, we set #3 (r) = 0, defining a point of 
isotropic conductivity. In Figure 5 we represent the fiber 
orientation (r) with respect to the local coordinate axes 
given by («(r),?i(r),? 2 W). 

If we have a uterine volume such that C is parallel to the 
z axis, then (r) can be written as a function of n(r) as 
follows: 



# 3 (r) = (P cos (a) 
where 



■ Fsin (a))n(r), 





~a 0 0 




"0 0" 


p = 


0 a 0 


, F = 


-HO 




0 0 -I/* 




0 0 0_ 


a = - 


V/(r) 


— , and 



(26) 



(27) 



V(VJ(r)) 2 + (V/(r)) 2 
|V/(r)|| 



(V^(r)) 2 + (V/(r)) 2 



(28) 



with Vy the ;-th component of the gradient. In the case of 
a spherical myometrium, (r) is given as by 



a 3 (r) = 



23/ cos a 

< s /x 2 +y 2 
R 



+ 



y sin a 



2^x 2 +y 2 
xs'ma 

+Jx 2 +y 2 



(29) 




Figure 5 Simplified illustrations of az (r) with respect to the local 
coordinates axis given by {w(r),Ti (r) f %(r)}. The blue plane 
contains the vectors n{r),k(r), andTi (r), and it is perpendicular to 
the gray plane formed by vectors % (r),T| (r) and «3(r).The orange 
plane is the cross section of the uterus perpendicular to the vector 
^§p-\ .The gray curve is the curve of symmetry rc (t) with rc (U ) 
and rc (t 2 ) extreme points of the curve. 



where R = ^x 2 + y 2 + z 2 . Note that, for a = 0, the main 
axis of the fibers runs vertically from the fundus to the 
cervix. 



Estimating myocyte fiber conductivities To the best 
of our knowledge, values of the intracellular and extra- 
cellular conductivity tensors have not been reported for 
the human myocyte, and therefore, these have to be esti- 
mated. To model the extracellular conductivity values a e \ 
and <7 e t, we assume a grid-type distribution of myocytes 
in the myometrium and use an estimate of the extracel- 
lular conductivity the human myometrium obtained by 
applying Archies law [23]. We describe myocytes as long 
cylinders with diameter d ce \\ and axis length / ce n, such 
that d ce \\ <^ / ce n Assuming that myocytes are uniformly 
arranged in a cubical grid with length lj = l ce \\ + 2A e and 
whose cross section has sides dj = d ce \\ + 2A e , then we 
have cr e \ and cr et as follows: 



cr e i = <7 e 



Get 



\ 2 \ 



and 



1 _ d ce \\l ce \\ \ 
djlj ) 



(30) 



(31) 
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where cr e is the conductivity of the extracellular medium in 
the myometrium. cr e can be computed using the effective 
myometrium conductivity <7m> available in the literature, 
and Archies law [23] as follows: 



a-p) 



m ' 



(32) 



where p is the fraction of the volume occupied by the 
myocytes and collagenous fibers in the tissue, and m is the 
so-called cementation factor, which depends on the shape 
and orientation of the myocyte in the tissue. 

To compute the intracellular conductivities <T[\ and 
we assume that the intracellular and extracellular domains 
have equal anisotropy ratios, that is, 



(33) 



and thus we need to compute g. We obtain an ana- 
lytical expression for g using reported values of the 
propagation speed of a transmembrane potential wave- 
form traveling on isolated tissue strips from pregnant 
human myometrium at term [28]. In particular, replac- 
ing (33), (12), and (13) in the bidomain equations (10) 
and (11), and solving v m (r,t) for a traveling wave solu- 
tion v m (§ • r—c t), where £ is a unitary vector pointing 
along the main axis of the myocyte and c is the speed of 
propagation, we obtain the following expression for g: 



£ =< 



2 c 2 eia m ci 



a e i k (v* - 2 v* + vlY 



(34) 



where g (x) = j^. Further, v\, v|, and V3 are the roots of 
the following polynomial in v m : 



V ■ (g + 1) G;V0 e (r,r) = -V ■ G' e Vv m (r,t), in M. 

(38) 

The above simplification is also known as the mon- 
odomain approximation of the bidomain equations, 
which, under suitable boundary conditions, allows the 
computation of v m and, thus,/ s , independent from 0 e . 

To set up boundary conditions for computing elec- 
trical potentials, we need to take into account the vol- 
ume conductor geometry (see Figure 3). In particular, we 
have two bidomain-monodomain interfaces: one between 
the myometrium M and abdominal volume A and one 
between the myometrium and the intrauterine cavity 
IA. Therefore, we have the following boundary condi- 
tions: (i) continuity of the interstitial potential 0 e at 
the perimetrium surface dM to the abdomen potential 
(j) A , (ii) flow of the normal component of / that crosses 
over from the uterus to the abdominal medium, (iii) no 
flow of the normal component of / s to the abdominal 
medium, (iv) continuity of the interstitial potential 0 e at 
the endometrium surface dU to the intrauterine cavity 
potential (pu, (v) flow of the normal component of / that 
crosses over from the uterus to the intrauterine cavity 
filled with amniotic fluid, (vi) no flow of the normal com- 
ponent of / s to the intrauterine cavity, (vii) no flow of 
the normal component of / that crosses over from the 
abdominal cavity to air, and (viii) either no flow or else 
flow of the normal component of / that crosses over from 
the intrauterine cavity, filled with amniotic fluid, to the 
fetus, depending on if it is covered with vernix caseosa 
(X = 0) or not (X ^ 0) [36]. These boundary conditions 
are summarized as follows: 



/ Om) = (Vm - n) (V 2 V m ) (v m - V 3 )-— (/3v mr + 8) , 

ky 

(35) 

with v mr denoting the resting transmembrane potential of 
the human myocyte. Note that in order to have g > 0, e\ 
must satisfy the following inequality: 



<t>e(r,t) = (pA(r,t), in dM, 



(39) 



0 < 61 < 



a el k (vj - 2 v* + v*) 2 
2 c 2 ci m 



Monodomain approximation and boundary conditions 

Assuming an equal anisotropy ratio, Eq. (33), simplifies 
the solution of the bidomain equations (10) and (11) by 
decoupling them as follows: 



V • -G e Vv m (r,£) = a m [c m — — +/i on 



(<T + 1) 



/ dv m (r,t 

T m at 

-/stimO, t)j in M, 



n M ' (G(V0i(r, t) + G;V0 e (r, t)) = n M • G^V0^(r, t), 

in dM, (40) 
n M ' GjVv m (r, t) = 0, in dM, (41) 

(36) 4>e(r,t) = 4>u(r,t), in dU (42) 
n u - (G(V0i(r, t) + G;V0 e (r, t)) = n u • G u V<hi(r, t), 

in dU, (43) 
n u -G[Vv m (r,t) = 0, in dU, (44) 
n A • G A V(j) A (r, t) = 0, in 3 A (45) 
njr • G w V0 w (r, t) = X (n T • G^V0^(r, t)), 

in dT, (46) 

(37) where nj is the normal vector to the surface ; in each case. 
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Summary of modeling assumptions 

In this section we list all the assumptions that were con- 
sidered while developing our model: 

• Volume conductor geometry: it is given by four 
spherical compartments. The uterine wall is a single 
layerwith uniform thickness. The abdominal, 
intrauterine, and fetal compartments have isotropic 
and homogeneous conductivity. The uterine 
compartment is anisotropic and homogeneous. All 
boundaries between compartments, except for the 
external boundary of the abdominal and fetal 
compartments, are electrically conductive. The fetal 
compartment boundary can be set to be electrically 
insulated or electrically conductive to take into 
account the presence of vernix caseosa. 

• Myometrium: It is a continuous medium formed by 
array of myocytes with same transmembrane 
potential shape across the whole domain. The 
extracellular and intracellular tissue conductivities 
are anisotropic, homogeneous, and proportional to 
each other (equal anisotropic ratio). The fiber 
orientation with respect to the main symmetry axis is 
kept fix across the whole domain. 

• Myocyte: It is a cell of cylindrical shape. Its electrical 
conductivity is assumed to have a cylindrical 
symmetry with higher conductivity along the main 
axis of the cylinder than along the axis denning the 
cross section. Transmembrane potential shape is 
plateau-type, and any myocyte can play the roll of 
pacemaker or follower by means of a stimulus current. 

Numerical computation of v m (r, t), 0 (r, t), and B (r, t) 

The computations of v m (r, t), 0 (r, t), and B(r, t) are given 
by the following procedure: 

• Step 1: Solve for v m (r, t) using Eqs. (37), (12), (13), 
and (14) subject to boundary conditions (41) and 
(44), and to initial conditions given by 



V • G A V(j)(r,t) = 0, in A 
V-G£,V0(r,f) =0, inW, 



(52) 
(53) 



v m (r, 0) = v mr , 

w m (r,0) = ( -v mr + 
\Y 

3v m (r,0) n , 

= 0, and 

dt 

3w m (r, 0) 



dt 



: 0. 



(47) 
(48) 

(49) 
(50) 



• Step 2: Solve for 0 e (r, t) in M and 0 (r, t) in A and 
U, using the solution of v m (r, i), computed in Step 1, 
and the following expressions: 

V • ( ^ + 1) G;V0 e (r, t) = - V • G;Vv m (r, t), in M, 

(51) 



subject to the boundary conditions (39), (40), (42), 
(43), (45), and (46). 
• Step 3: Solve for B(r, t) using Eq.(l), and computing 
the total current density /(r, t) in the whole uterine 
domain, using the solutions of v m (r, f), 0 e O"> 0» an d 
0(r, t), obtained in Steps 1 and 2, and considering the 
following boundary conditions: 

n T x (Bjr(r f t) - B u (r, t)) = 0, in d (54) 

n u x (B M (r, t) - B u (r, t)) = 0, in dU, (55) 

n M x (B A (r, t) - B M (r, t)) = 0, in dM, (56) 

n A x (B A (r, t) - B e (r, t)) = 0, in dA, (57) 

n A xB s (r,t) =0,mdS, (58) 

where £ is an additional volume that must be 
incorporated in order to compute the magnetic field 
at the abdominal surface using FEM. The boundary 
conditions at the interfaces 3 J 7 , dU, dM, and dA 
describe the continuity of the magnetic field between 
domains, and the boundary condition at 3£ , the 
external boundary of £ , is of electric isolation. 

To compute the solution in each of the above steps, we 
use the FEM solver COMSOL Multiphysics, running on a 
server with eight 64-bit processors at 2.3 GHz, with 32 GB 
RAM. 

Results 

In the following, we illustrate our modeling approach by 
considering the electrophysiological characteristics of the 
myometrium and a spherical volume conductor geometry 
In Figure 6, we illustrate a spherical, four-compartment 
volume conductor geometry used in our numerical exam- 
ple. We define a spherical myometrium with a 16 cm 
radius measured from the center to dM. and assume the 
uterine wall has a uniform thickness of 1 cm. We also con- 
sider a spherical fetus with a 12 cm radius concentric to 
the myometrium and fully covered with vernix caseosa, 
i.e., X = 0 in (46). The abdominal compartment is also 
spherical, with a 21 cm radius shifted —3 cm from the cen- 
ter of the myometrium in the x axis. This shift is to take 
into account that the uterus is closer to the abdominal sur- 
face than to the dorsal surface. We set the coordinate axis 
of reference at the center of the myometrium. 

The conductivity values for each compartment are given 
in Table 1. In particular, to compute the extracellu- 
lar myometrial conductivity tensors, we use the average 
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values for the uterine myocyte dimensions at term based 
on data reported in [1,4,9] (see Table 2). We assume the 
average human myocyte shape to be a long cylinder with a 
small cross section; therefore, we use a cementation factor 
m = 4/3 (see [23] for more details on the computation of 
this factor). The volume fraction p occupied by myocytes 
and collagenous fibers in the myometrium is set to 0.6. 
In order to consider an average myometrial fiber architec- 
ture that contains both circularly- and obliquely-oriented 
fibers we choose the fiber orientation angle a to be 45°. 
Figure 7 illustrates the global structure of the myometrial 
fiber orientation for this angle. 

We select the model parameters of the ionic current 
model using phase-space analysis, using as a reference the 
average plateau-type transmembrane potential recorded 
from isolated tissue strips of human myometrium at term 
[13,30]. In particular, the average resting potential, con- 
sidering the results reported from the 37th week of preg- 
nancy onwards, is approximately —56 mV. The plateau has 
an average depolarization of — 27zbl mV that terminates in 
0.9 ±0.2 minutes by an abrupt repolarization to the resting 



level [13,30]. Table 3 gives the parameter values used in 
the numerical example. Note that we compute the surface- 
area to volume-ratio a m using the myocyte dimensions in 
Table 2. 

In [19], the authors reported that, in the human uterus, 
there may be a preferential direction of propagation of 
contractions, and thus of transmembrane potential prop- 
agation, from the fundus toward the isthmus, which 
could aid in the expulsion of the fetus. Therefore, in 
order to study this assumption with our model, we con- 
sider / st i m with N p = 1, vi = 2 V, and h\{r,t) = 
{1, if 0 < t < 0.1 , 0.15 < ||r|| < 0.16 and z > 0.15 ; 0, 
otherwise. The size and intensity of the pacemaker area 
are chosen in order to obtain a stable traveling wave- 
form solution to the bidomain equations on the spherical 
myometrium. 

To solve the set of equations that model the myome- 
trial current-source densities and their electromagnetic 
fields at the level of the abdomen, we use a three-step 
procedure along with Finite Element Methods (FEM) 
(See the Methods Section for more details). The FEM 



Table 1 Conductivity values of the volume conductor Table 2 Myocyte dimensions and Archied law parameters 



geometry 






Symbol 


Value 


Reference 


Symbol 


Value 


Reference 


dcell 


7 jLtm 


[1] 


Ga. 


0.2 S/m 


[36] 


/cell 


450 /xm 


[1] 


Gu 


1 .74 S/m 


[36] 


dj 


8 jam 




Gt 


0.2 S/m 


[36] 


lj 


451 jLtm 




d e | 


0.68 S/m 


Eq. (30) 




0.5 S/m 


[36] 


Oet 


0.22 S/m 


Eq. (31) 


m 


4/3 


[4,23] 


5" 


0.8 


Eq. (34) 


P 


0.6 


[5,23] 
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-0.2 -0.2 

Figure 7 Geometry and fiber orientation in spherical myometrium given by a = 45°. 



discretization of the whole volume conductor is done 
using tetrahedral elements. The elements used to dis- 
cretize each domain consists of 75, 964 elements divided 
as follows: 1, 650 elements in the fetus, 7, 605 ele- 
ments in the intrauterine cavity, 8,299 elements in the 
myometrium, and 19, 036 elements in the abdominal cav- 
ity. To compute the magnetic field using FEM, an addi- 
tional domain (domain £ in (58)) concentric to the volume 

Table 3 Ionic current model parameters 



Symbol Value Reference 

c m 0.01 F/m 2 [36] 

v mr -0.056 V [30] 

a m 5.75871 0 5 m~ 1 Table 2 

61 200 co m 2 

6 2 0.09 1/s 
V] -0.02 V 
v 2 -0.04 V 
v 3 -0.065 V 
k 1 0 4 1 A/ 2 
8 0.0520 V 
y 0.1 

P 1 

c 1.15 cm/sec [28] 



conductor geometry has to be added. Specifically, we add a 
concentric sphere with radius of lm. The discretization of 
this domain consists of 39, 374 elements. The discretized 
boundaries consists of 10,996 triangular elements, and 
the number of degrees of freedom of the whole domain 
is 159, 130. In all the Steps of the numerical computation 
procedure described above we use the generalized min- 
imal residual method (GMRES) for solving the resulting 
system of linear equations after FEM discretization. The 
backward differentiation formula is used to discretize the 
time derivative in Step 1. 

Figure 8 shows several snapshots of the FEM solu- 
tion for one pacemaker on the fundus of a spheri- 
cal myometrium, assuming anisotropy as illustrated in 
Figure 7. Figure 8 (a) -(c) illustrates the transmembrane 
potential and source current density distribution at the 
myometrium, Figure 8 (d)-(f) shows the electrical poten- 
tial at the abdominal surface, and Figure 8 (g)-(i) illustrates 
the magnetic field density at the abdominal surface. The 
magnetic field measured at the abdominal surface, Bmmg 
(the magnetomyogram field), is proportional to B nA , the 
projection of B onto the normal vector of the abdominal 
surface, n^. 

In Figure 9 (a) we illustrate the temporal response of 
the FEM solutions for the transmembrane potential at 
different elevations over time. It can be seen that a sta- 
ble traveling waveform has been established as the shape 



La Rosa etal. BMC Medical Physics 2012, 12:4 
http://www.biomedcentral.eom/1 756-6649/1 2/4 



Page 13 of 16 



Time = 10 [sec] 




(a) 



Time = 36 [sec] 



y ftxit En) 



(b) 



Time — £E [sec] 



0.15 
u 1 

* 0 

;: 

m 

-u 1 
0. 15 



r 



Time = 10 fseel 



Time = 10 fsecl 



-0.1 0 C-.1 




y - i x x t [n] 

(c) 



E 

I 



C 2 
J. LB 
C 1 

H 

7 

-c.i 
-j is 

-C 5 





(d) 

Time = 36 [sac] 



(g) 

Time = 36 [sac] 




Figure 8 FEM solution at time instants t = 1 0 [s], 36 [s], 55 [s] for one pacemaker on the fundus of a spherical myometrium, assuming 
anisotropy. (a)-(c) transmembrane potential and source current density distribution at the myometrium, (d)-(f) electrical potential at the 
abdominal surface, and (g)-(i) magnetic field density at the abdominal surface. 



remains the same. Also, the maximum depolarization 
is —16 mV, the average potential in the plateau area 
is —25 mV, and the transmembrane potential duration, 
before hyperpolarization, is approximately 35 s, which are 
fair approximations with respect to the average recorded 
transmembrane potentials discussed in [13,30]. Note that 
our ionic current model introduces hyperpolarization, 
which constrains the excitability of the cell, and thus con- 
secutive contractions can only take place until v m reaches 
resting potential. In our case, our model can simulate a 
minimum period of one contraction every 240 s, allowing 



us to model labor contractions whose period reduces to 
about one contraction every 300s. 

In Figure 9(b) we illustrate the percentage of the con- 
tracting myometrial volume as function of time, which in 
[7,9] was used as a reference to compute the changes in 
the intrauterine pressure due to a contraction. Interest- 
ingly, we observe that the percentage of the contracting 
myometrial volume as a function of time has the sym- 
metric properties and time duration of the intrauterine 
pressure waveforms of a pregnant human myometrium at 
term, as discussed in [9]. 
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100 



Time [aec] Time aec 



(a) (b) 

Figure 9 (a) Temporal response of the FEM solutions for transmembrane potential at different elevations; (b) Percentage of contracting 
myometrial volume as a function of time. 



Discussion 

Our numerical simulation results (see Figure 8) show that, 
because of the anisotropy in the conductivity, the direc- 
tion of the current density / s is rotated at a certain angle 
from the main direction of the propagating transmem- 
brane potential v m , and it is the transversal component 
of this current, parallel to the x-y plane, which generates 
the magnetic field B nA . This observation is in agreement 
with the analysis presented in [37], and it is important to 
take it into account when interpreting the magnetic field 
measurements generated by uterine contractions in the 
presence of volume conductor geometry. Therefore, the 
spatial signature of B nA is highly dependent on the fiber 
orientation of the myometrium. Because of the proxim- 
ity of the sensors and the myometrium, it is not strictly 
applicable to assume a moving dipole parallel to the direc- 
tion of propagation of the transmembrane potential as the 
main model for the current source generating the mea- 
sured magnetic field. This last interpretation might be 
suitable in the case where the transverse length of the 
transmembrane potential front is short in comparison to 
the area covered by the array of sensors. In contrast, if 
the transverse length is larger and thus not covered by 
the measured area, (for example, when several cells are 
recruited), then it is suitable to consider a moving line 
source (stretched ring) model instead. 

Though our volume conductor geometry is an oversim- 
plification of a real anatomical structure, certain aspects 
of it stand as a fair approximation of the volume conduc- 
tor geometry when performing MMG measurements. For 
example, the SARA system, known as the superconduct- 
ing quantum interference device array for reproductive 



assessment a (SQUID Array for Reproductive Assessment) 
is a passive, stationary, floor-mounted instrument at 
which the patient sits and leans her abdomen against a 
concave surface which contains an array of 151 sensors. 
The effect of leaning the abdomen on the concave sur- 
face works as a way to standardize the abdominal surface 
making the spherical model very suitable to represent 
the measuring surface. Also, a spherical uterus is good 
approximation at the early stage of pregnancy but not nec- 
essarily through all stages. In this sense, a more realistic 
population-average uterine model should be constructed 
from magnetic resonance images. Unfortunately, expos- 
ing of pregnant patients to MRI is not currently the norm, 
and thus having access to average uterine geometry for 
different stages of pregnancy is not yet possible. 

The ionic current model based on extended FHN 
equations can reproduce a good approximation of 
the plateau-type transmembrane potential recorded in 
human myocytes, however, it introduces hyperpolariza- 
tion, which constrains the excitability of the cell, and 
thus consecutive contractions can only take place until 
v m reaches resting potential. In our case, our model can 
simulate a minimum period of one contraction every 240 
s, allowing us to model labor contractions whose period 
reduces to about one contraction every 300 s. In addition, 
note that a larger 62 value can extend the transmem- 
brane potential duration to values closer to the average 
duration reported in [13,30]. However, a larger 62 value 
also extends the duration of hyperpolarization and there- 
fore the plateau of the curve describing the percentage of 
contracting myometrial volume. This last observation is 
the result of the combined effect of the transmembrane 
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potential duration and the geometry of the uterus. Clin- 
ical observations on the shape of the uterine pressure 
wave have found correlation between the rising slope of 
the waveform and contraction efficiency [38], namely, the 
steeper the slope the more efficient the contraction is. 
Additional numerical examples, using one pacemaker at 
the same position, show that our model is consistent with 
the latter observation, since we found that changing the 
propagation speed of the activity modifies the shape of 
the percentage of the contracting myometrial volume as a 
function of time. In particular, a faster speed of propaga- 
tion of the main transmembrane potential front increases 
the rising slope the percentage of the contracting myome- 
trial volume, however, above a certain threshold it reduces 
the duration of the contraction. Placing two pacemakers 
at each extreme of the spherical domain, one at the fundus 
and the other one at the cervix, doubles the rising slope 
of the percentage of the contracting myometrial volume 
curve, and, depending on the duration of the transmem- 
brane potential, it introduces a plateau in the curve, which 
implies that all volume is contracting. The symmetry of 
the contracting myometrial volume curve generated using 
our model is primarily due to the symmetry of the volume 
conductor used in our simulations and the duration of the 
transmembrane potential simulated. 

Modeling the different stages of pregnancies can be 
accommodated using this model. First the volume con- 
ductor geometry should be scaled to the average vol- 
ume size corresponding to stage of pregnancy of inter- 
est. Specifically, the fetal, the uterine, and the abdominal 
shapes should be modified accordingly using reported 
values in the literature (e.g., see [3] for uterine shapes). 
Second, the conductivity values of the volume conductor 
model should be set with respect to the stage of preg- 
nancy. For example, the conductivity of the intrauterine 
cavity given by the amniotic fluid is known to change 
as the pregnancy develops (e.g., see [36]). Regarding the 
equal anisotropy ratio, this can be computed by using the 
reported values of the speed of propagation, the resting 
transmembrane potential, and the updated ionic current 
parameters such these fit the transmembrane potential 
for the specific stage. Third, the boundary conditions 
remain the same across all periods, noting that for bound- 
ary condition given by Eq. (46), X should be set to 1 or 
0 depending on the presence of vernix caseosa. Fourth, 
incorporation of a bursting type transmembrane potential 
can be done using the FHN model by adding a periodic 
function to Eq. (13), which control the excitability of the 
cell. However, more realistic and complex shapes of trans- 
membrane potential can be designed by incorporating, for 
example, the ionic current model presented in [17]. 

Additional approaches to validate our mathematical 
model should consider comparisons against real mag- 
netic and potential field measurements at the level of the 



abdomen, which leads automatically to the next natural 
question: how to solve the inverse problem of our model? 
Explicitly, given a set of measurements and a paramet- 
ric model, can we infer the value of certain parameters of 
interest, such as the number of stimui and their positions, 
the conductivity tensor in the domain, initial and bound- 
ary conditions, etc? Solving the inverse problem stands as 
the connecting bridge between multiscale modeling of the 
pregnant uterus and clinical applications. In particular, as 
a first approach we are planning to use MMG measure- 
ments to estimate the current density in the myometrium 
and thus search for features that characterize contractions 
patterns which lead to preterm-labor. 

Conclusion 

We proposed a multiscale-forward electromagnetic 
model of uterine contractions during pregnancy. Our 
model incorporates knowledge of the electrophysiologi- 
cal aspects of uterine contractions during pregnancy at 
both the cellular and organ levels. We applied a bido- 
main approach for modeling the propagation of the 
myometrium transmembrane potential v m on the uterus 
and used this to compute the action potential 0 and the 
magnetic field B at the abdominal surface. We introduced 
a modified version of the FitzHugh-Nagumo equation for 
modeling the ionic currents in each cell. Though our ionic 
current model does not consider Ca 2+ dynamics explic- 
itly, the simplicity of the modified equation allows for the 
propagating action potential to be modeled under well 
denned conditions as shown in this paper. We also pro- 
posed a general approach to design conductivity tensors 
in the myometrium and to estimate the conductivity ten- 
sor values in the extracellular and intracellular domains. 
We introduced a simplified geometry for the problem 
and proposed a discretized model solution based on a 
finite element method approach. Finally, we illustrated 
our modeling approach through a numerical example 
by modeling uterine contractions at term. Our model is 
potentially important as a tool for helping characterize 
contractions and for predicting labor using MMG and 
EMG. 

As part of our future work, we will investigate pear- 
shaped uterine domains as a way to approximate better 
the uterine geometry. We will also include more accurate 
ionic current models as in [13-15,17] and will consider 
spatial variations of the fiber orientation. 

Endnote 

a SARA was built in collaboration with VSMMedTech Ltd., 
Canada and is installed at the University of Arkansas for 
Medical Sciences (UAMS) Hospital. 
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